/* SPDX-License-Identifier: SunMicrosystems */
/* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. */

/**
 *
 * This family of functions implements the cubic root of :math:`x`.
 *
 * Synopsis
 * ========
 *
 * .. code-block:: c
 *
 *     #include <math.h>
 *     float cbrtf(float x);
 *     double cbrt(double x);
 *     long double cbrtl(long double x);
 *
 * Description
 * ===========
 *
 * ``cbrt`` computes the cubic root of the input value.
 *
 * Mathematical Function
 * =====================
 *
 * .. math::
 *
 *    cbrt(x) \approx \sqrt[3]{x}
 *
 * Returns
 * =======
 *
 * ``cbrt`` returns the cubic root of the input value.
 *
 * Exceptions
 * ==========
 *
 * Does not raise exceptions.
 *
 * Output map
 * ==========
 *
 * +---------------------+--------------+---------------------+--------------+--------------+---------------------+--------------+--------------+
 * | **x**               | :math:`-Inf` | :math:`<0`          | :math:`-0`   | :math:`+0`   | :math:`>0`          | :math:`+Inf` | :math:`NaN`  |
 * +=====================+==============+=====================+==============+==============+=====================+==============+==============+
 * | **cbrt(x)**         | :math:`-Inf` | :math:`\sqrt[3]{x}` | :math:`x`                   | :math:`\sqrt[3]{x}` | :math:`+Inf` | :math:`qNaN` |
 * +---------------------+--------------+---------------------+--------------+--------------+---------------------+--------------+--------------+
 *
 *///

#include <math.h>
#include "../common/tools.h"

#ifndef __LIBMCS_DOUBLE_IS_32BITS

static const uint32_t
B1 = 715094163, /* B1 = (682-0.03306235651)*2**20 */
B2 = 696219795; /* B2 = (664-0.03306235651)*2**20 */

static const double
C  =  5.42857142857142815906e-01, /* 19/35     = 0x3FE15F15, 0xF15F15F1 */
D  = -7.05306122448979611050e-01, /* -864/1225 = 0xBFE691DE, 0x2532C834 */
E  =  1.41428571428571436819e+00, /* 99/70     = 0x3FF6A0EA, 0x0EA0EA0F */
F  =  1.60714285714285720630e+00, /* 45/28     = 0x3FF9B6DB, 0x6DB6DB6E */
G  =  3.57142857142857150787e-01; /* 5/14      = 0x3FD6DB6D, 0xB6DB6DB7 */

double cbrt(double x)
{
#ifdef __LIBMCS_FPU_DAZ
    x *= __volatile_one;
#endif /* defined(__LIBMCS_FPU_DAZ) */

    int32_t hx;
    double r, s, t = 0.0, w;
    uint32_t sign;
    uint32_t high, low;

    GET_HIGH_WORD(hx, x);
    sign = hx & 0x80000000U;     /* sign= sign(x) */
    hx  ^= sign;

    if (hx >= 0x7ff00000) {
        return (x + x);    /* cbrt(NaN,INF) is itself */
    }

    GET_LOW_WORD(low, x);

    if ((hx | low) == 0) {
        return (x);    /* cbrt(0) is itself */
    }

    SET_HIGH_WORD(x, hx);   /* x <- |x| */

    /* rough cbrt to 5 bits */
    if (hx < 0x00100000) {    /* subnormal number */
        SET_HIGH_WORD(t, 0x43500000);   /* set t= 2**54 */
        t *= x;
        GET_HIGH_WORD(high, t);
        SET_HIGH_WORD(t, high / 3 + B2);
    } else {
        SET_HIGH_WORD(t, hx / 3 + B1);
    }


    /* new cbrt to 23 bits, may be implemented in single precision */
    r = t * t / x;
    s = C + r * t;
    t *= G + F / (s + E + D / s);

    /* chopped to 20 bits and make it larger than cbrt(x) */
    GET_HIGH_WORD(high, t);
    INSERT_WORDS(t, high + 0x00000001, 0);


    /* one step newton iteration to 53 bits with error less than 0.667 ulps */
    s = t * t;    /* t*t is exact */
    r = x / s;
    w = t + t;
    r = (r - t) / (w + r); /* r-s is exact */
    t = t + t * r;

    /* retore the sign bit */
    GET_HIGH_WORD(high, t);
    SET_HIGH_WORD(t, high | sign);
    return (t);
}

#ifdef __LIBMCS_LONG_DOUBLE_IS_64BITS

long double cbrtl(long double x)
{
    return (long double) cbrt((double) x);
}

#endif /* #ifdef __LIBMCS_LONG_DOUBLE_IS_64BITS */
#endif /* #ifndef __LIBMCS_DOUBLE_IS_32BITS */
